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Abstract 

A hybrid Maxwell solver for fully relativistic and electromagnetic (EM) particle- 
in-cell (PIC) codes is described. In this solver, the EM helds are solved in 
k space by performing an FFT in one direction, while using hnite differ¬ 
ence operators in the other direction(s). This solver eliminates the numerical 
Cerenkov radiation for particles moving in the preferred direction. Moreover, 
the numerical Cerenkov instability (NCI) induced by the relativistically drift¬ 
ing plasma and beam can be eliminated using this hybrid solver by applying 
strategies that are similar to those recently developed for pure FFT solvers. 

A current correction is applied for the charge conserving current deposit 
to correctly account for the EM calculation in hybrid Yee-FFT solver. A 
theoretical analysis of the dispersion properties in vacuum and in a drifting 
plasma for the hybrid solver is presented, and compared with PIC simulations 
with good agreement obtained. This hybrid solver is applied to both 2D and 
3D Cartesian and quasi-3D (in which the helds and current are decomposed 
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into azimuthal harmonics) geometries. Illustrative results for laser wakefield 
accelerator simulation in a Lorentz boosted frame using the hybrid solver in 
the 2D Cartesian geometry are presented, and compared against results from 
2D UPIC-EMMA simulation which uses a pure spectral Maxwell solver, and 
from OSIRIS 2D lab frame simulation using the standard Yee solver. Very 
good agreement is obtained which demonstrates the feasibility of using the 
hybrid solver for high hdelity simulation of relativistically drifting plasma 
with no evidence of the numerical Cerenkov instability. 

Keywords: PIC simulation, hybrid Maxwell solver, relativistic plasma 
drift, numerical Cerenkov instability, quasi-3D algorithm 


1. Introduction 

Fully relativistic, electromagnetic particle-in-cell (PIC) codes are widely 
used to study a variety of plasma physics problems. In many cases the solver 
for Maxwell’s equations in PIC codes use the hnite-difference-time-domain 
(FDTD) approach where the corresponding differential operators are local. 
This locality leads to advantages in parallel scalability and ease in implement¬ 
ing boundary conditions. However, when using PIC codes to model physics 
problems, including plasma based acceleration [1] in the Lorentz boosted 
frame, relativistic collisionless shocks [ai, and fast ignition parti¬ 

cles or plasmas stream across the grid with speeds approaching the speed of 
light. In these scenarios, the second order FDTD Maxwell solvers support 
light waves with phase velocities less than the speed of light. This prop¬ 
erty of the FDTD solver leads to numerical Cerenkov radiation from a single 
particle that is moving near the speed of light. In addition, when beams or 
plasmas are moving near the speed of light across the grid a violent numerical 
instability known as the numerical Cerenkov instability (NCI) arises due to 
the unphysical coupling of electromagnetic modes and the Langmuir modes 
(main and higher order aliased beam resonance). The beam resonances are 
at a; -|- 27r/r/At = {ki + 27ri/i/Axi)r’o, where /i and ui refer to the time and 
space aliases. At and Axi are the time step and grid size, and the plasma is 
drifting relativistically at a speed Vq in the 1-direction. 

The NCI was hrst studied more than 40 years ago [7]. However, it has re¬ 
ceived much renewed attention [HlElliniElinillalE] since the identification 
[I5l[l6j of this numerical instability as the limiting factor for carrying out 
relativistic collisionless shock simulations [SIE], and Lorentz boosted frame 
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simulations [IZl HSl HH |2U] of laser wakefield acceleration (LWFA) [T] . 

This early and recent work on the NCI [3 El El CHI HIl HSl El] have 
shown that the NCI inevitably arises in EM-PIC simulations when a plasma 
(neutral or non-neutral) drifts across a simulation grid with a speed near 
the speed of light. Analysis shows that it is due to the unphysical coupling 
of electromagnetic (EM) modes and Langmuir modes (including those due 
to aliasing). As a result, significant recent effort has been devoted to the 
investigation and elimination of the NCI so that high fidelity relativistic 
plasma drift simulations can be routinely performed Uni El ini ESI EH Eoj. 

In previous work Ref. IHED], we examined the NCI properties for the 
second order Yee solver [22] , as well as a spectral solver [231 [23] (in which 
Maxwell’s equation are solved in multi-dimensional k space). We note that 
what we refer to as simply a spectral solver, others refer to as a pseudo- 
spectral time domain (PSTD) solver. The NCI theory developed in [H EH] 
were general and it could be applied to any Maxwell solver. It was found 
that in the simulation parameter space of interest the fastest growing NCI 
modes of these two solvers are the (/r, z/i) = (0, ±1) modes, where [x and ui 
dehned above are the temporal aliasing, and spatial aliasing in the drifting 
direction of the plasma. The (/r, z/i) = (0, ±1) modes for both solvers reside 
near the edge of the fundamental Brillouin zone (for square or cubic cells), 
and can be eliminated by applying a low-pass hlter. However, due to the 
subluminal EM dispersion along the direction of the drifting plasma in the 
Yee solver, the main NCI mode (/r, z/i) = (0,0) of the Yee solver has a 
growth rate that is of the same order as its (/i, z^i) = (0,±1) counterpart, 
and these modes reside close to the modes of physical interest. However, 
the (/i, z^i) = (0,0) NCI mode in the spectral solver has a growth rate one 
order of magnitude smaller than the (/x, z^i) = (0, ±1) modes due to the 
superluminal dispersion of spectral (EFT based) solver. Furthermore, as 
shown in [10] these (/r, ui) = (0, 0) modes can be moved farther away from the 
physics modes and their harmonics by reducing the time step in the spectral 
solver, and can be fully eliminated by slightly modifying the EM dispersion 
in the spectral solver. Using these methods, it was demonstrated in ra that 
a spectral EM-PIC can perform high fidelity simulations of relativistically 
drifting plasmas where the LWFA physics is highly nonlinear with no evidence 
of the NCI. 

A multi-dimensional spectral Maxwell solver has a superluminal disper¬ 
sion relation in all the propagation directions. This is due to the fact that 
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the first order spatial derivatives in the Maxwell’s equation are greater than 
iV-th order accurate (where N is the number of grids) since we are solving the 
Maxwell’s equation in k space. This superluminal dispersion relation leads 
to highly localized (p., = (0,0) NCI modes and the reduction of their 

growth rates (compared with their Yee solver counterpart). In this paper, 
we propose to use a hybrid Yee-FFT solver, in which the FFT is performed 
in only one direction, namely the drifting direction of the plasma, while keep¬ 
ing the hnite difference form of the Yee solver in the directions transverse to 
the drifting direction. In other words, EM waves moving in the 1 direction 
will have a superluminal dispersion (due to the N-th order accurate spatial 
derivatives) while those moving in the 2 (and 3 in 3D) directions will have 
a subluminal dispersion due to the second-order-accurate spatial derivatives. 
The advantages of this approach over a full FFT solver is that the held solver 
is local in the transverse directions so that better parallel scalability than a 
fully FFT based solver can be achieved (assuming the same parallel FFT 
routines are used). In addition, it is easier to include a single FFT into 
the structure of mature codes such as OSIRIS |26]. Furthermore, this idea 
works well with a quasi-3D algorithm that is PIC in r — z and gridless in 0 
[271128] . where the FFTs cannot be applied in the f direction. We note that 
recently a method for achieving improved scalability for FFT based solvers 
was proposed [23] in which FFTs are used within each local domain, but it 
introduces as yet unquantihed errors in the longitudinal helds. The relative 
advantages and tradeoffs between the variety of approaches being proposed 
will be better understood as they begin to be used on real physics problems. 

We use the theoretical framework for the NCI developed in Ref. [9l[in| to 
study the NCI of the hybrid solver. As we show below, the fastest growing 
NCI modes for the proposed hybrid solver behave similarly to those for the 
spectral solver. In k space they reside at the edge of the fundamental Bril- 
louin zone for square or cubic cells. More importantly, the (/r, z/i) = (0, 0) 
NCI mode for the hybrid solver has almost the same properties (pattern, 
growth rates) as that of a spectral solver. The NCI can therefore be effi¬ 
ciently eliminated in the hybrid solver by applying the same strategy as in 
the spectral solver. Moreover, simulations have shown that the NCI proper¬ 
ties of the quasi-3D r — z PIC and gridless in 0 algorithm [IZIEH] are similar 
to that of 2D Cartesian geometry [29|. Therefore, the idea of a hybrid Yee- 
FFT solver can be readily applied to quasi-3D geometry. We also note that 
the use of local FFTs in domains along z [25] could be also be used within 
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the hybrid approach described here. 

In this paper, we hrst discuss the algorithm for the hybrid Yee-FFT 
Maxwell solver in section In section we apply the theoretical framework 
in Ref. [9l [10] to study the NCI properties of the hybrid solver analytically 
and in PIC simulations. We compare OSIRIS [2n| results with the hybrid 
solver against UPIC-EMMA [TT] results with a fully spectral (FFT based) 
solver. In section]^ it is shown that the strategies used to eliminate the NCI 
for purely spectral solvers are also valid for the hybrid solver. In section 
we extend the hybrid solver idea to the quasi-3D algorithm in OSIRIS and 
present simulation studies of the NCI properties in this geometry. We then 
present 2D OSIRIS simulations of LWFA in a Lorentz boosted frame using 
the new hybrid solver. Very good agreement is found when comparing simu¬ 
lation results using the hybrid solver in OSIRIS against results from 2D lab 
frame OSIRIS using Yee solver and 2D Lorentz boosted frame UPIC-EMMA 
[TT] simulations using spectral solver. Last, in section]^ we summarize the 
results and mention directions for future work. 

2. Hybrid Yee-FFT solver 

The basic idea of the hybrid Yee-FFT solver is that the theoretical frame¬ 
work developed in [91 do] indicates that the NCI is easier to eliminate when 
EM waves are super luminal along the direction of the plasma drift. This can 
be accomplished with higher order solvers or with an FFT based solver in 
the drifting direction of the plasma (denoted as 1-direction). We note that it 
is more difficult to satisfy strict charge conservation (Gauss’s law) for higher 
order hnite difference solvers. Here we replace the hnite difference operator 
of the hrst spatial derivative d/dxi in the Maxwell’s equation in Yee solver 
with its FFT counterpart that has an accuracy greater than order N. We 
then correct for this change in the current deposit to maintain strict charge 
conservation. Without loss of generality, in the following we will briehy de¬ 
scribe the algorithm of the Yee-FFT solver in two-dimensional (2D) Cartesian 
coordinate. The straightforward extension to the 3D Cartesian case is also 
discussed. 

2.1. Algorithm 

We start from the standard algorithm for a 2D Yee solver, in which the 
electromagnetic helds E and B are advanced by solving Faraday’s Law and 
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where the EM held E and B, and cnrrent j are dehned with the proper half¬ 
grid offsets according to the Yee mesh izg. If we perform a Fonrier transform 
of Eq. Q-(|^ in both Xi and X 2 , and in time, Maxwell’s eqnations rednce to 
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In vacuum where j = 0, the corresponding numerical dispersion relation for 
the EM waves is 


+ [k]i) 


( 10 ) 


The idea of a hybrid Yee-FFT solver is to keep the hnite difference opera¬ 
tor [k ]2 = sin(A: 2 Aa; 2 / 2 )/(Ax 2 / 2 ) in the directions transverse to the drifting 
direction, while replacing the hnite difference operator [k]i in the drifting di¬ 
rection with its spectral counterpart [k]i = ki. To achieve this, in the hybrid 
solver we will solve Maxwell’s equations in ki space. The current is deposited 
locally using a rigorous charge conserving scheme that is equivalent to [30] . 
For the EM held and current, we hrst perform an FFT along xi so that all 
helds are dehned in {ki,X 2 ) space. After that we apply a correction to the 
current in the drifting direction 
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where ji is the corrected current. In [25], the current is also corrected where 
they combine a pure FFT solver with a charge conserving current deposit. 
This correction ensures that Gauss’s Law is satished throughout the duration 
of the simulation if it is satished initially, as will be discussed in more detail 
in section 12.31 After the current correction we advance the EM held as 
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where ki = 2 ttki/N and N is the number of grids in xi direction, and 
Ki = 0,1,..., N/2 — 1 is the mode number. Note in the hybrid solver, the 
EM helds E, B, and current j have the same temporal and spatial centering 
as in the Yee solver, and 


= exp 



(18) 


is the phase shifting due to the half grid offsets of the Ei, i? 2 , 3 , and ji in the 
1-direction. Compared with the standard Yee solver algorithm, it is evident 
that if we replace —iki with the corresponding hnite difference form we can 
recover the standard 2D Yee algorithm. 


2.2. Courant condition 

The Courant condition of the hybrid solver can be easily derived from 
the corresponding numerical EM dispersion Eq. (10). Substituting into Eq. 
(10) the hnite difference operator in time [ca] 
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Note the k range of the fundamental Brillouin zone is \ki\ < tt/Axi, \k 2 \ < 
TT/Ax 2 , we can obtain the Courant limit on the hybrid solver 
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For square cells with Axi = Ax 2 , this reduces to At < 0.537Axi. 











2.3. Charge conservation 

In the hybrid Yee-FFT solver, we rely on the Faraday’s Law and Am¬ 
pere’s Law to advance the EM held. On the other hand, the local charge 
conserving current deposition Isa ensures the second-order-accurate hnite 
difference representation of the continuity equation. 
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where is an arbitrary scalar quantity. Therefore, when combining this 
scheme with the second order accurate Yee solver. Gauss’s Law is rigorously 
satisfied at every time step if it is satisfied at t = 0. However, when the 
hybrid solver is used together with the charge conserving current deposition 
scheme, we need to apply a correction to the current, as shown in Eq. 0. 
in order that the Gauss’s Law is satisfied at every time step. This can be 
seen by first performing Fourier transform in the Xi direction for Eq. ([2^, 
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then applying the divergence operator of the hybrid solver to the left and 
right hand side of the Ampere’s Law, Eq. (15)-(17). Using Eq. (25), we 
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is satisfied at t = 0, it is satisfied at each time step. 
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2.4- 3D Cartesian geometry 

It is straightforward to extend the hybrid solver to 3D cartesian geometry. 
In 3D Cartesian coordinates, we solve Maxwell’s equation in {ki, X 2 ,x^) space 
where we use the same second order accurate hnite difference form of the Yee 
solver in the 2 and 3 directions. As in the 2D Cartesian case, the current 
correction is applied to ji to ensure the Gauss’s Law is satished. We have 
implemented the hybrid solver in 2D and 3D with current correction in our 
hnite-difference-time-domain (FDTD) code OSIRIS [2S]- 


3. Numerical Cerenkov instability 


To investigate the NCI properties of the hybrid solver, we hrst consider 
its corresponding numerical dispersion relation. Employing the general the¬ 
oretical framework established in Ref. 13 EQ], we can calculate in detail the 
NCI modes for any Maxwell solver. The roots of the numerical dispersion 
relation that lead to the NCI can be found numerically by solving Eq. (17) 
in |9], or by the analytical expression in Eq. (19) of [10]. For convenience we 
present the corresponding numerical dispersion and analytical expressions of 


Eq. (17) of [9] in Appendix A For the Yee solver the k space representation 
of the hnite difference operator is 


[k]i 


sin(/cjAa;j/2) 

^Xi/2 


(28) 


where i = 1, 2 in 2D. Meanwhile, in the hybrid solver the k space operator in 
the drifting direction is replaced with that of the spectral solver [fc]i —)• ki. 
By substituting the respective operators for each direction into Eq. (19) of 
Ref. HD] [or Eq. 


(A.4) in Appendix A 
In Fig. [l 


we can rapidly hnd the set of NCI 
a)-(d), we plot the (/x, z^i) = (0,0) 


modes for the hybrid solver 
and (/i, z^i) = (0, ±1) modes for the hybrid and spectral solvers by scanning 
over the {ki,k 2 ) space in the fundamental Brillouin zone and solve for the 
growth rates of the corresponding unstable modes. The parameters used to 
generate this plot are listed in Table 

We can see from Fig. (a) and (b) that the (/i, ui) = (0, ±1) NCI 
modes of the two solvers reside near the edge of the fundamental Brillouin 
zone, although the patterns are slightly different due to their different hnite 
difference operators in the ^direction, which leads to the slightly different 
EM dispersion curves. In Fig. (e) and (f) we show how different EM 
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Parameters 

Values 

grid size {kpAxi, kpAX 2 ) 

( 0 . 2 , 0 . 2 ) 

time step UpAt 

OAAxi 

boundary condition 

Periodic 

simulation box size {kpLi, kpL 2 ) 

51.2x51.2 

plasma drifting Lorentz factor 

7 = 50.0 

plasma density 

n/np = 100.0 


Table 1: Crucial simulation parameters for the 2D relativistic plasma drift simulation. Up 
is the reference plasma density, and = iirq^Up/me, kp = ujp (c is normalized to 1). 

dispersion curves leads to different (/t, ui) = (0, ±1) NCI modes for the two 
solvers. These modes are distinct, and far removed from the modes of physical 
interest, and are relatively easy to eliminate. 

More importantly, we see from Fig. (c) and (d) that the hybrid solver 
leads to (/i, = (0,0) NCI modes that are very similar to their spectral 

solver counterpart. The pattern of the (/i, ui) = (0, 0) modes for these two 
solvers are both four dots (in 2D) and highly localized in the fundamental 
Brillouin zone. We also use the theory to perform parameter scan to study the 
dependence of growth rates (of the fastest growing mode) and the locations 
in k space of the (//, zz^) = (0,0) modes on At/Axi for the hybrid solver, 
and compare this result against that of the fully spectral solver, as shown 
in Fig. (a) and (b). We likewise carried out OSIRIS simulations using 
the hybrid solver and UPIC-EMMA m [31] using the spectral solver, to 
compare against theoretical results. Very good agreement is found between 
theory and simulations. Fig. (a) and (b) show that both the ki location, 
and growth rates of the (/r, ui) = (0, 0) modes are almost identical for the two 
solvers. This indicates that, just like the spectral solver, the growth rate of 
the (yU, ui) = (0,0) modes of the hybrid solver is reduced, while their location 
in ki increases when the time step is reduced. 

In Fig. 1^ (c) and (e) we show the locations of the unstable (/i, ui) = 
(0,±1), and (/i, zzi) = (0,0) NCI modes for the hybrid solver in OSIRIS for 
2D Cartesian geometry. The agreement between Fig. (c) and Fig. 0(b). 
and between Fig. (e) and Fig. (d) are excellent. 

The main advantage of the purely spectral solver regarding its NCI prop¬ 
erties in comparison to a purely FDTD solver is that the superluminal dis¬ 
persion relation makes it much easier to eliminate the NCI modes at (/i, ui) = 
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(0, 0): the modes have a growth rate that is one order of magnitude smaller 
than that for the (/i, vi) = (0, ±1) modes, their locations are highly localized 
in k space, and they can be moved away from the modes of physical interest 
by reducing the time step. We showed above that similar NCI properties can 
be achieved by using a hybrid FDTD-spectral solver, where the Maxwell’s 
equation are solved in Fourier space only in the direction of the plasma drift. 
Comparing with an EM-PIC code using a multi-dimensional spectral solver 
which solves Maxwell’s equation in k space, there are advantages when solv¬ 
ing it in {ki,X 2 ) space in 2D [and {ki,X 2 ,X 3 ) space in 3D]. Firstly, the hybrid 
solver saves the FFT in the other directions; secondly, since the solver is 
FDTD in the directions transverse to the drifting direction, it is easier to 
integrate the algorithm into existing FDTD codes such as OSIRIS where the 
parallelizations and boundary conditions in the transverse direction can re¬ 
main untouched. Last but perhaps most important, the idea that one can 
obtain preferable NCI properties by solving Maxwell’s equation in ki space 
in the drifting direction can be readily extended to the quasi-3D algorithm 
as we can solve the Maxwell’s equation in (fci,p,-0) space. 

4. Elimination of the NCI modes 

In Ref. [To], we proposed strategies to eliminate the NCI in the spectral 
solver. These strategies can be readily applied to the hybrid solver. For 
square (or cubic) cell, the pattern of the fastest growing modes resides in a 
narrow range of ki near the edge of the fundamental Brillouin zone. Therefore 
we can apply a low-pass hlter in ki to the current to eliminate the fastest 
growing modes. Since the fields are already in ki space when solving the 
Maxwell’s equations, the hltering can be done efficiently by applying a form 
factor to the current only in ki. 

As for the (/i, ui) = (0, 0) mode, if they are near the main or higher or¬ 
der harmonics of the physical modes, we can move them away and reduce 
their growth rates by simply reducing the time step. To further mitigate the 
{fi, Ui) = (0, 0) NCI modes when they are far away from the physical modes, 
one can modify the EM dispersion relation, according to the procedure de¬ 
scribed in Ref. [10], to completely eliminate them. In Fig. we plot how 
the modification is accomplished in the hybrid solver. As shown in Fig. 

(a) except for the bump region for most ki the [fcji for a particular ki is ki 
itself; near the bump, the [fcji for ki is ki + Akmod, where Ak^od is a function 
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of ki with 


k\ kim TT 
kll kirn 2 

where ku, ki^ are the lower and upper ki to be modihed, kim = {kii + kiu)/2, 
and Akmod,ma.x is the maximum value of Ak^od- The values of ku, and 
Akmod,raax are determined by the position of the (/x, z^i) = (0,0) modes and 
their growth rates. According to the NCI theory, for the parameters in Table 
1^ when the [k]i is as dehned in Fig. g (a) (with kn/kgi = 0.15, klu/kgl 
0.26, and Akmod,max/kgi = 0.01), there is no unstable (/i, ui) = (0,0) NCI 
modes, i.e., the (/i, ui) = (0,0) mode has a theoretical growth rate of zero. To 
verify the theoretical results in the hybrid solver, in Fig. g(b) we plot the E 2 
energy growth with and without the modihcation. In these simulations we 
used the parameters in Table g The blue curve in Fig. g (b) represents the 
case without the modihcation, while the red and black curves are those with 
the modihcation to ki. The cases with blue and red curves used quadratic 
particle shapes, while the case for the black curve used cubic particle shapes. 
We have likewise plotted the E 2 spectra at the time point t = 3200 oJq^ 
indicated in Fig. g (c) and (d) for the two cases with the modihcations 
(red and black curves in g (b)). We can see from Fig. g (b) and (c) that 
after the modihcation, the growth rate of the (/x, vi) = (0, 0) NCI modes 
reduces to zero. Meanwhile, the red curve rises later in time due to the 
(/X, i^i) = (±1,±2) NCI modes. As we showed in Ref. [101 the growth rate 
of these higher order modes can be reduced by using higher order particle 
shape. Therefore when cubic particle shapes are used, as is the case for the 
black curve, the (/x,z/i) = (±1, ±2) NCI modes do not grow exponentially 
and are therefore much less observable in the corresponding spectrum at 
t = 3200 in Fig. g(d) as compared to|(c). 

5. hybrid solver in quasi-3D algorithm 

As mentioned in section g the idea of the hybrid solver can be easily 
incorporated into the quasi-3D algorithm [271 EH] in which the helds and 
current are expanded into azimuthal Fourier modes. We can obtain the 
hybrid Yee-FFT solver for the quasi-3D algorithm by using FFTs in the z (xi) 
direction and hnite difference operators in r (^ 2 ) direction in the equations for 
each azimuthal mode. Note in quasi-3D OSIRIS we use a charge conserving 


Akmod Akmod,raa.x COS 
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current deposition scheme for the Yee solver (as described in [2B]), therefore 
for the hybrid solver adapted for the quasi-3D algorithm we can apply the 
same current correction for the use of FFTs to ji in order that the Gauss’s 
Law is satished throughout the duration of the simulation. 

The NCI properties of the hybrid solver for the quasi-3D algorithm are 
similar to that of the 2D Cartesian geometry [25]. While a rigorous NCI 
theory for the quasi-3D algorithm is still under development, we can empir¬ 
ically investigate the NCI for this geometry through simulation. In Fig. 

(d) and (f) we plot the Er data at a time during the exponential growth 
of the EM helds due to the NCI, which shows the (/i, z^i) = (0, ±1) and 
(p, vi) = (0, 0) modes for the hybrid solver in quasi-3D geometry. For the 
Er data, we conduct an FFT in xi and a Hankel transform in X 2 - Similarly 
to the 2D Cartesian case, we isolate the (/r, z/i) = (0, 0) modes by applying 
a low-pass hlter in the current in ki space to eliminate the fastest grow¬ 
ing (/i, z^i) = (0, ±1) NCI modes. The parameters used in the simulations 
are listed in Table and a conducting boundary is used for the upper r 
boundary. We kept azimuthal modes of m = —1, 0,1 in the simulations. 

By comparing Figs. (c)-(f) it can be seen that the pattern of the NCI 
modes are similar for the {x 2 , a;i) and (r, z) geometries. We have also plotted 
the dependence of the growth rate and ki position of the (/i, ui) = (0,0) NCI 
modes for the quasi-3D geometry in Fig. (a) and (b). These plots show 
that when the time step decreases the growth rates of the (p, z/i) = (0, 0) NCI 
modes in the quasi-3D geometry decreases, while the ki position increases 
(and move away from the physical modes), in a nearly similar fashion to 2D 
Cartesian geometry. This indicates that the same strategies for eliminating 
NCI in 2D Cartesian geometry can be applied to the quasi-3D geometry. 
The fastest growing modes residing at the edge of the fundamental Brillouin 
zone can be eliminated by applying a low-pass hlter in the current. The 
(/i, vi) = (0, 0) NCI modes can be mitigated by either reducing the time step 
to lower the growth rate and move the modes away from the physics in ki 
space, or by modifying the [k]i operator as discussed in section]^ to create 
a bump in the EM dispersion along the ki direction. We have implemented 
the modihcation to the [k]i operator into the hybrid solver for the quasi-3D 
OSIRIS code, and have conhrmed that this modihcation completely eliminate 
the = (0,0) NCI modes. The coefficients used for the modihcation 

are the same as those for the 2D Cartesian case discussed in section HI 
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6 . Sample simulations 

In this section, we present preliminary results of Lorentz boosted frame 
LWFA simulations using the hybrid solver in OSIRIS. For comparison, we 
performed simulations with the same parameters using UPIC-EMMA which 
uses a spectral Maxwell solver. Table lists the simulation parameters. We 
use a moving antenna in both cases to launch lasers into the plasma. The 
results are summarized in Fig. 

In Fig. 1^ (a)-(b) the Ei held at t' = 3955a;(^^ for simulations with both 
the hybrid solver and spectral solver in the Lorentz boosted frame are plotted, 
where ujq is the laser frequency in the lab frame. Both the spectral solver and 
hybrid solver give similar boosted frame results, and there is no evidence of 
NCI affecting the physics in either case. We plot the line out of the on-axis 
wakeheld in Fig. (c), which shows very good agreement with one another. 
The very good agreement can also be seen when we transformed the boosted 
frame data back to the lab frame. In Fig. (d)-(f) we plot the on-axis Ei 
held for the OSIRIS lab frame data, the transformed data for the OSIRIS 
boosted frame simulation with the hybrid solver, and the transformed data 
from theUPIC-EMMA boosted simulation at several values of time in the lab 
frame. As seen in Fig. (d)-(f), the transformed data from the two boosted 
frame simulations agrees very well with each other. Note the displacement 
of the lineouts between the lab frame data and boosted frame data is due to 
the different group velocity of the laser between the Yee, spectral, and hybrid 
solver. Since the hnite difference operator [k]i is the same for the spectral and 
hybrid solver, the group velocity of the laser along its propagation direction 
is the same for the two boosted frame simulations. As a result the on-axis 
transformed laser data of the two boosted frame simulations almost reside 
on top of one another and are more accurate. 

7. Summary 

We proposed to use a hybrid Yee-FFT and a rigorous charge conserving 
current deposit for solving Maxwell’s equations in order to eliminate the nu¬ 
merical Cerenkov instability in PIC codes when modeling plasmas or beams 
that drift with relativistic speeds in a particular direction. In this solver 
we solve the Maxwell’s equation in ki space along the drifting direction {xi 
direction), and use second order hnite difference representation for the deriva¬ 
tives in the other directions. This provides greater than Y-th order accuracy 
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Plasma 

density uq 

1.148 X 10 ^no7fe 

length L 

7.07 X lO^kQ^/jb 

Laser 

pulse length r 

70Mko^jb{l + /3b) 

pulse waist W 

117.81ko^ 

polarization 

3-direction 

normalized vector potential Oq 

4.0 

2D boosted frame simulation 

grid size Axi ^2 

Om82ko^jb{l + /3b) 

time step At/Axi 

0.225 

number of grid {'jb = 14) 

8192x512 

particle shape 

quadratic 


Table 2: Parameters for a 2D LWFA simulations in a Lorentz boosted frame that were 
used for in 2D Cartesian geometry with the hybrid solver in OSIRIS and with a fully 
spectral solver in UPIC-EMMA. The laser frequency Wq ^nd number fco in th® inb frame 
are used to normalize simulation parameters. The density is normalized to the critical 
density in the lab frame, no = TOeWo/(47re^). 


for the spatial derivatives in the Xi direction, while keeping the locality of 
the held solve and current deposit in the directions transverse to 1. For the 
current deposit, we start from the charge conserving deposit in OSIRIS and 
then correct it so that it still satishes the continuity equation for the hybrid 
solver. Thus, Gauss’s law remains rigorously satished at every time step if it 
is satished initially. 

It is found from the NCI theory that such a hybrid solver has similar 
NCI properties in comparison to a full spectral solver that solves Maxwell 
equation in multi-dimensional k space. As a result, the (p., = (0,0) NCI 

modes have a growth rate one order of magnitude smaller than the fastest 
growing (/x, z/i) = (0, ±1) NCI modes, and are highly localized. In addition, 
the growth rates of the (/i, z/i) = (0, 0) modes decrease as one reduces the 
simulation time step, and their locations in Fourier space also move farther 
away from the physics. 

Compared with the spectral solver, the hybrid solver performs an FFT 
only along the drifting direction of the plasma. As a result, it saves the 
computation of FFT in the other directions if this ultimately becomes an 
issue for parallel scalability. In addition, it can be readily adapted into fully 
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operational FDTD codes without the need to modify various boundary condi¬ 
tions in the transverse directions. Very importantly, this idea can be readily 
applied to the quasi-3D algorithm in which the quantities are decomposed 
into azimuthal harmonics. In this algorithm FFTs cannot be used in the f 
direction. We demonstrate the feasibility of the hybrid Yee-FFT solver in 
2D/3D Cartesian geometry, as well as in the quasi-3D geometry. Although 
we have not conducted a rigorous theoretical analysis for the NCI in the r-z or 
quasi-3D geometries, we hnd in simulations that the hybrid solver in quasi- 
3D geometry has very similar NCI properties to that in the 2D Cartesian 
geometry. 

We show that the strategy to eliminate NCI in the hybrid solver for 
2D/3D Cartesian geometry, as well as quasi-3D geometry, is similar to that 
for the spectral solver. The fastest growing NCI modes can be eliminated 
by applying a low-pass hlter in the current. The (/i, z/i) = (0, 0) NCI modes 
can be eliminated by reducing the time step which both reduces their growth 
rates and moves them away from the physical modes in Fourier space. These 
NCI modes can also be fully eliminated by slightly modifying the the EM 
dispersion relation along ki direction at the location in Fourier space where 
the (/i, vi) = (0, 0) modes reside. This approach is demonstrated in both 
Cartesian and quasi-3D geometry. 

We showed that the new hybrid solver in OSIRIS can be used to conduct 
2D LWFA simulations in a Lorentz boosted frame. With the low-pass hlter 
applied to current and using reduced time step, we observe no evidence of NCI 
affecting the physics in the simulation. Very good agreement is found between 
the results from OSIRIS with the hybrid solver, UPIC-EMMA simulations, 
as well as OSIRIS lab frame simulations with the standard Yee solver. This 
demonstrates the feasibility of using the hybrid solver to perform high hdelity 
relativistic plasma drift simulation. 

This work was supported by US DOE under grants de-sc0008491, de- 
SC0008316, DE-FC02-04ER54789, DE-FG02-92ER40727, by the US National 
Science Foundation under the grant ACI 1339893, and by NSFC Grant 
11175102, thousand young talents program, and by the European Research 
Gouncil (ERG-2010-AdG Grant 267841), and by LLNL’s Lawrence Fellow¬ 
ship. Simulations were carried out on the UGLA IIoffman2 and Dawson2 
Glusters, and on Hopper cluster of the National Energy Research Scientihc 
Gomputing Genter. 
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Appendix A. Numerical dispersion for relativistically drifting plasma 
and NCI analytical expression in hybrid solver 


According to Ref. [a [in], the numerical dispersion for the hybrid solver 
can be expressed as 


[ojr-[k]Ei[k]m-[k]E2[k]B2-^{-ir 

7 


Sj2iSE2[(^] — 5's3[A:]_Bino) \ 


u' — k[vo 


+ C = 0 

where C is a coupling term in the dispersion relation 
(-1)" 


(A.l) 


C = 


7 


SjiSElUj'[k]E2[k]B2iVo — 1) + Sj2SE2[k]E2[k]B2i^' — k'^Vo) 

+ Sji[k]E2{SE2[k]Bik2Vo ~ *S'B3A;2no[i^])l (^-2) 


and for the hybrid solver 


[k]El — [k]Bl — ki [k]E2 — [k]B2 — 


sin(/c2Aa;2/2) 

Aa;2/2 


(A.3) 


We can expand u' around the beam resonance u' = in Eq- (A.l), 
and write u' = k'^v^ + 5u', where 5u' is a small term. This leads to a cubic 
equation for 5oj' (see [101 for the detailed derivation), 

^2^0;'^ + + C25uj' + T>2 = 0 (A.4) 

where 


A2 =2eo'6 

B 2 ='Co|'Co ~ [k]El[k]Bl — [k]E2[k]B2 - ■^{ — B)^Sj2{SE2il — ^*553 [fcjei) | 

^^2 =‘^(-l)^|eo5',2(CoA^3Wi?l - SE2io) - ilS,Ak]E2k2SE2[k]Bl 
+ io[k]E2{Sj2SE2[k]B2 — A;2Cl*S'53.^o) | 

D2 ="^{-mo[k]E2k2S^, (^SE2[k]Bi - (A.5) 
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where 


= cos(fciAt/2) 

Cl = — sin(fciAt/2)At/2 


Co 

Co 

ki 


sin(/ciAt/2) 
At/2 

cos(fciAt/2) 

kl + Vlkgl - 


We use 


/ sin(/c^Aa:j/2) y^^ 
V Axi/2 ) 


(A.6) 


(A.7) 


as well as use the corresponding interpolation functions for the EM helds 
used to push the particles 


SeI — S/,lSZ,2Si,3( —1)^^ Se2 — Sz,lS/,2SZ,3 

SbI = S/,iSz,2Sz,3 Sb2 = Sz,lSi^2Sz,3( —1)^^ 


SeI — SZ,iSz,2Sz,3 

SbZ = ■5z,iSz,2Sz,3(~1)^^ 

(A.8) 


when using the momentum conserving field interpolation, and use 

SeI = ■SZ-l,lSi^2Sz,3( — Se2 = Sz,iSZ-1,2Sz,3 Sei = S;^iS;^2'SZ-1,3 

SbI = Si^iSi-i^2Sl-l,3 Sb2 = Sz_i4Si,2SZ-l,3( —1)^^ Sb3 = SZ-l,lSZ-l,2Sz,3( —1)*^^ 

(A.9) 

when using the energy conserving held interpolation. The (—1)^^ term is due 
to the half-grid offsets of these quantities in the 1 direction. With respect to 
the current interpolation, 


Sjl — SZ-l,lSZ,2SZ,3( —1)^^ Sj2 — Sz,lSi_i,2SZ,3 Sj3 — Sz,lSZ,2SZ-l,3 


(A.IO) 


We note that we use expressions for charge conserving current deposition 
scheme that are strictly true in the limit of vanishing time step At —)■ 0. The 
coefficients A 2 to D 2 are real, and completely determined by ki and k 2 . By 


solving Eq. (A.4) one can rapidly scan the NCI modes for a particular set of 
(/x,z/i). 
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spectral solver hybrid solver 

Instability mode (mu,nut) = 0,1 Instability mode (mu,nu1) = 0,1 



k1/kg1 k1/kg1 

Instability mode (mu,nu1) = 0. 0 Instability mode (mu,nu1) = 0, 0 





Figure 1: The pattern of the (/x, vi) = (0, ±1) modes for the two solvers are shown in (a) 
and (b). The pattern of the (Hji'i) = (0,0) modes for two solvers are shown in (c) and 
(d). The intersection between the EM dispersion relations with the first spatial aliasing 
beam modes for the full spectral solver and the hybrid solver are shown in (e) and (f). 
When generating these plots we use Axi = Ax 2 = 0.2 and At = 0.08 ■ Other 

parameters are listed in Table 
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Figure 2: In (a) and (b) the dependence of the growth rate and ki for the fastest growing 
(/r, ui) = (0,0) mode on the time step is shown. The four lines correspond to the theoretical 
prediction for the hybrid solver in 2D, results from OSIRIS and UPIC-EMMA simulations 
for the spectral and hybrid solvers in 2D Cartesian geometry, and results for the hybrid 
solver in the quasi-3D geometry (where the k 2 is obtained from a Hankel transform). In 
(c)-(f) the spectrum of E 2 {Ep) is plotted for OSIRIS simulations with the hybrid solver in 
2D Cartesian or the quasi-3D geometry. In (c) and (d) results from runs where no filter in 
ki is used to eliminate the (/i, vi) = (0, ±1) modes. In (e) and (f) a filter in ki was used to 
eliminate the (/i, vi) = (0, ±1) modes and now the (^, vi) = (0, 0) modes are seen. These 
results show that the 2D Cartesian and quasi-3D geometries have very similar properties 
and that the strategies used to eliminate the NCI in 2D Cartesian can be applied to the 
quasi-3D case. 
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modification of [k]i operator 



E: energy evoiution 



quadratic particie shape [red curve in (a)] 


cubic particie shape [biack curve in (a)] 



k, /k . 

1 gi 


k, / k , 

1 gi 


Figure 3: In (a) the perturbation to [A:]i that is used to eliminate the (/r, ui) = (0,0) NCI 
modes is shown. In (b) the evolution of the logj^Q for a reference case and for two 
cases with the EM dispersion modification (one with quadratic and another with cubic 
particle shapes). In (c) and (d), the spectrum of E 2 at t = 3200 is shown for the 
two cases with the EM dispersion modifications. In (c) quadratic particle shapes are used, 
while in (d) cubic particle shapes are used. 
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Figure 4: Comparison between OSIRIS lab frame, OSIRIS with the hybrid solver in the 
boosted frame and UPIC-EMMA in the boosted frame. In (a) and (b), 2D plots of Ei 
for OSIRIS with the hybrid solver and UPIC-EMMA at t' = 3955a;(^^ are shown in the 
boosted frame, where ujq is the laser frequency in the lab frame. In (c), lineouts along the 
laser propagation direction of the same data are shown. In (d)-(f), lineouts of the Ei data 
transformed back to the lab frame are shown. The colored lines correspond to an OSIRIS 
lab frame simulation, an OSIRIS hybrid solver simulation in the Lorentz boosted frame, 
and UPIC-EMMA simulation in the Lorentz boosted frame. 
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